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Abstract: We present a first application of a previously published method for the com- 
putation of QCD processes that is accurate at next-to-leading order, and that can be 
interfaced consistently to standard shower Monte Carlo programs. We have considered Z 
pair production in hadron-hadron collisions, a process whose complexity is sufficient to test 
the general applicability of the method. We have interfaced our result to the HERWIG 
and PYTHIA shower Monte Carlo programs. Previous work on next-to-leading order cor- 
rections in a shower Monte Carlo (the MC@NLO program) may involve the generation of 
events with negative weights, that are avoided with the present method. We have compared 
our results with those obtained with MC@NLO, and found remarkable consistency. Our 
method can also be used as a standalone, alternative implementation of QCD corrections, 
with the advantage of positivity, improved convergence, and next-to-leading logarithmic 
accuracy in the region of small transverse momentum of the radiated parton. 
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1. Introduction 

Next-to-leading order (NLO) QCD calculations have by now become the standard for the 
study of processes relevant to collider physics. It is well known that the inclusion of NLO 
terms are important for a reliable estimate of cross sections. One finds corrections of 
the order of 30 to 100% for typical production processes. Until recently, the results of 
NLO computations were not included in shower Monte Carlo (SMC) models, since their 
inclusion is highly non trivial. In ref. [1] a method (referred to as MC@NLO) was proposed 
to include the results of NLO calculations in a parton shower Monte Carlo, and it was 
applied to several processes^ [2, 3] in connection with the HERWIG implementation [4]. 
The MC@NLO method in general must be adapted to the SMC algorithm being used. 

^For the full list of processes implemented in MC@NLO, see the MC@NLO web page 
http: //www.hep.phy . cam. ac .uk/theory/webber/MCatNLO/. 
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Furthermore, the event generation is not strictly positive, i.e. negative weighted events are 
generated. 

In ref. [5] a novel method was proposed, that overcomes these problems, and allows 
one to compute the generation of the hardest radiation (i.e. the radiation with the largest 
transverse momentum) independently of the particular SMC employed for the subsequent 
shower. It was found there that, in the framework of angular ordered SMC's, in order 
to preserve the correctness of the soft radiation spectrum, it is necessary to include the 
soft radiation coherently emitted from the partons produced in the hardest splitting. Such 
coherent radiation was named "truncated shower" , since it is a shower that stops when the 
angular variable approaches the angle of the hardest emission. 

In ref. [5], only the general strategy was discussed, while the path to follow for a 
practical implementation was not studied. The generation of the hardest event was only 
sketched there, and no indication was given on how to actually perform it in practice. 
Furthermore, the implementation of the truncated shower in a SMC program was left to 
future work. This last problem should not however be overemphasized. It is relevant only 
to angular ordered SMC's that fully enforce soft coherence. At present, only the HERWIG 
program does this. Other SMC's implementations do not care to enforce coherence at the 
same level of accuracy.^ Furthermore, it may be possible to perform the generation of the 
truncated shower independently of the particular SMC used. 

The problem of interfacing NLO calculations to SMC's can thus be split into three 
independent issues: 

1. Generation of the hardest emission: one constructs an algorithm to generate the 
hardest emission with NLO accuracy. The resulting events can be put in a standard 
form, such as the "Les Houches Interface for User Processes" (LHI) [8]. 

2. Generation of the truncated shower: given the hardest event, on the LHI, one can add 
the truncated shower and put the event back in the LHI. This step is only required 
if one wants to maintain the correct soft emission pattern, and is using an angular 
ordered SMC. 

3. Showering, hadronization, decays: any SMC that complies with the LHI requirements 
can now perform the rest of the showering, provided it also implements a p-j- veto. 
Today's most popular SMC's satisfy these requirements. 

Following the above strategy, the problem becomes more manageable, since it no longer 
requires to modify or rewrite a full SMC implementation. Furthermore, several independent 
solutions to each of the above points may be pursued by independent researchers. In the 
present work, we tackle item 1 of the above list, in the specific case of Z boson pair 
production in hadronic collisions. We have chosen this process for the following reasons: 

• It is an important process for LHC physics. 

^There arc programs that generate soft radiation using dipole type formalisms [6, 7], but they do not 
treat initial state radiation consistently. 
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• It involves initial state radiation, which is more difficult than final state radiation. 



• It is readily extended to the similar WZ and WW production processes. 

• It should be easily generalizable to other important processes, like, for instance, single 
boson, Higgs and heavy flavour production. 

We will now recall the basic formulae for the generation of the hardest emission given 
in ref. [5], and illustrate in brief the most relevant features of the method that we have 
developed. 

In ref. [5], it was required that the phase space kinematics is factorized in terms of 
Born (v) and radiation (r) variables. In the case of ZZ production the Born variables can 
be chosen as the pair invariant mass Mzz and rapidity Yzz, and the cosine of a polar angle 
9i. In the case of two-body subprocesses, 6i is the angle formed by the direction of one 
of the outgoing Z's and the beam axis, in the rest frame of the ZZ system. In the case 
of 0{as) real-emission corrections, the ZZ system has non-zero transverse momentum. 
In this case we define a three-dimensional frame in the ZZ rest system, and define 9i as 
the angle of one of the Z^s relative to its third axis. There is some arbitrariness in the 
choice of this frame. The only important requirement is that in the limit of zero transverse 
momentum of the ZZ pair, its third axis should become parallel to the collision axis. A 
specific choice is described in detail in ref. [9]. The radiation variables (r) are chosen to 
be X = MJ^/s, where s is the invariant mass of the incoming partonic system, a variable 
— 1 < y < 1, equal to the cosine of the scattering angle of the emitted parton in the rest 
frame of the incoming parton system, and an azimuthal variable 02 which parametrizes the 
relative azimuthal angle of the ZZ system with respect to the radiated parton. With these 
definitions, when the r variables approach the soft (x ^ 1) or coUinear (y ±1) limits, 
the kinematics of the ZZ system approaches the corresponding Born kinematics with the 
given V variables. We denote by 

B{v)d^y, V{v)d^^, R{v,r)d<Pyd^r, C{v,r)d<Pyd^r (1-1) 

the Born, soft- virtual, real and counterterm contributions to the cross section, respectively. 
The differential cross section for the hardest emission can be written schematically as 



da = B{v)d^v 



A{v,0) + A{v,k^{v,r))^^^d^r 

B(v) 



(1.2) 



where 



B{v) = B{v) + V{v) + j d^r [R{v, r) - C{v, r)] (1.3) 



A{v,Pt) = exp 



- / ^^^9{k^{v,r)-p^)d^r 



(1.4) 



and k'Y{v,r) is the transverse momentum of the emitted parton. As written above the 
calculation of the probability of the hardest emission would seem computationally intensive. 



-3- 



In particular, the computation of B{v) requires one three-dimensional integral for each v. 
In order to circumvent this problem, we introduce the function 

B{v, r) = N[B{v) + V{v)] + R{v, r) - C{v, r) (1.5) 

where 
so that 

B{v) = / B{v,r)d<!>r . (1.7) 



Standard proccdiires are available to generate unweighted v, r events from the distribu- 
tion B{v,r)d^r d^v ■ This is exactly what wc need: we generate v,r values in this way, 
and we ignore the r value, which amounts to integrating over the r variables. Thus, the 
generation of unweighted events distributed according to B(v) is no more expensive, from 
a computational point of view, than generating unweighted events for the real emission 
matrix elements. The generation of the radiation variables r also looks computationally 
intensive, but it can be performed in an efficient way by using the veto method. 

In the following, we will illustrate all the details of the implementation of the procedure 
outlined above. In Section 2 we collect the kinematics and cross section formulae for the 
ZZ production process. In Section 3 we write down in full detail eq. (1.2). In Section 4 
we discuss important issues having to do with the choice of scales, and the accuracy in 
the Sudakov form factors for the generation of the hardest event. In Section 5 wc first 
describe how one would perform a straightforward implementation of unweighted event 
generation using the given hard cross section. In Subsections 5.1 and 5.2 we illustrate in 
detail our method for the generation of the Born variables v and of the radiation variables 
r. Explanations of the Monte Carlo techniques that we have used are reported in detail in 
the appendices, for ease of reference. 



2. Kinematics and cross section 

The differential cross section for the production of Z boson pairs in hadronic collisions was 
computed in ref. [9] up to order Ug. The effects of spin correlations of the decay products 
of the two vector bosons have been computed in refs. [10, 11], but will not be included here 
for simplicity. In this section, we formulate the result of ref. [9] in a form which is suitable 
for the generation of the hardest event using the procedure proposed in ref. [5]. 

The order-Q:^ cross section for the process H1H2 —>■ ZZ + X can be written as the sum 
of four terms: 

da = da^""^ + da^^^^ + da^^^ + da^"^ . (2.1) 

Here da^^^ is the leading-order (Born) cross section. The term da^^''^ collects order-ccg 
contributions with the same two-body kinematics as the Born term, namely one-loop cor- 
rections and real-emission contributions in the soft limit. Finally, dcr'-^^ represents the 
cross section for real emission in a generic configuration, while da^"^^ is a remnant of the 
subtraction of collinear singularities, and describes real emission in the coUinear limit. 
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At leading order, the only relevant parton subprocess is 

qipi) + q{p2)^ Zih)+Z{k2), (2.2) 

where g is a quark or antiquark of any flavour, and q the corresponding antiparticle. Particle 
four-momenta are displayed in brackets; we have Pi=P2 = 0, fe^ = fc^ = m|, where ruz is 
the Z boson mass. We introduce the usual Mandelstam invariants 

s = {pi+P2f, t = {pi-kif, U={pi-k2f, (2.3) 

related hy s + t + u = 2m|. 

Event generation is conveniently performed in terms of the invariant mass Mzz and 
the rapidity Yzz of the Z boson pair in the laboratory frame. They are given by 

Ml^ = {ki + k2f = xiX2S (2.4) 

V {Vl+V2f + {Pl+P2f 1, 

2 {pi + P2r - {Pi + P2r 2 X2 

where S is the squared center-of-mass total energy of the colliding hadrons, xi and X2 are 
the fractions of longitudinal momenta carried by the incoming partons, and we have used 
the fact that the ZZ pair has zero rapidity in the center-of-mass frame of the colliding 
partons. It follows that 



xi 



Af 2 / 1 

-^e^^^=xbi; X2 = ^ ^ e-^^^ = Xb2; dxidx2 = ^dYzzdMzz- (2.6) 



We adopt as two-body kinematic variables the set v = {Mzz, Yzz, cos 6i} (which we will 
call the Born variables henceforth), where 9i is the angle between pi and ki in the partonic 
center-of-mass frame, so that 

t = ml-^il-(3cos9,), (2.7) 

where 

Using eq. (2.6) and the usual expression of the two-body phase space measure d^2, it is 
immediate to check that 

d^2 dxi dx2 = — ^ d cos 9i dM^^ dYzz ■ (2.9) 
In order to keep our notation similar to that of ref. [5], we define 

d% = dcos Oi dMl^ dYzz- (2.10) 
The appropriate integration region for the variables v is 

1 1 
4m|<ML<5, -log-|^<y^^<--log-|^, -l<cos0i<l. (2.11) 
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The Born cross section is given by 



(2.12) 



where 



^?(^'/^) = (^62,/i)H?(ML,t). 



IGttS 



(2.13) 



The index q runs over all quarks and antiquarks; f^{x, jj) denotes the distribution function 
of parton q in the hadron iJ, and /x is a factorization scale. The function M.^qq{s,t) is 
the squared invariant amplitude, summed over final-state polarizations and averaged over 
initial-state polarizations and colors, divided by the relevant flux factor: 



2s Nc 



t u Amis A / I 1 

- + 7 + - "^M 72 + ^ 

u t tu \t'^ 



(2.14) 



where gqy and gq^A denote the vector and axial-vector couplings of the quark q to the Z 
boson, and Nc = 3 is the number of colours. 

Order-ccg contributions to the cross section arise from one-loop corrections to the two- 
body process eq. (2.2), and from rcal-cmission subprocesses at tree level. The contribution 
of one-loop diagrams must be summed to the one-gluon emission cross section in the soft 
limit, in order to obtain an infrared-finite result. The resulting contribution has the same 
kinematic structure as the leading-order term: 



(2.15) 



where 



(2.16) 



log- 



Ml. 



11^ 



(6 + 16 log P) + 32 log^ /3 - -TT^ 



and Cp = 4/3. The explicit expression of Ai^^^^^^\s, t, ji) is given in Appendix B of ref. [9]. 

The dependence in A^^^'"^*^^-* is implicit through its dependence upon a^. Notice also that 
the explicit dependence in eq. (2.16) is due to the collinear subtraction. For simplicity, 
we have chosen the same value for the renormalization and factorization scales. Since the 
Born process is of order in a^, there is no explicit ji dependence due to renormalization, 
and therefore the renormalization scale only appears as the argument of a^. 
Next, we consider the real-emission subprocesses 



q{pi) + Q{pi) 
q{pi) + g{p2) 
g{pi) + q{p2) 



Z{h) + Z{k2) + g{k) 
Z{ki)+Z{k2) + q{k) 
Z{ki)+Z{k2) + qik) 



(2.17) 
(2.18) 
(2.19) 
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in a generic kinematical configuration. The processes (2.17-2.19) are characterized by five 
independent scalar quantities, which we choose to be 

s = (Pl+P2)^ ife = Uk = {p2-kf, qi = {pi-kif, ^2 = (p2-A;2)^, (2.20) 

as in ref. [9]. We introduce the variables 

x = —^; y = cose, (2.21) 

where 9 is the scattering angle of the emitted parton in the partonic center-of-mass system. 
With these definitions, 

tk = ~{l - x){l - y); Uk = ~{l-x){l + y). (2.22) 

It is easy to show that in the case of the subprocesses (2.17-2.19) one has 

1 Xis + Uk 1 Xi2- (l-x)(l+y) 2 c /ooQ\ 
Yzz = -log — — = -log— - — -; Mzz = xxiX2S, (2.23) 

2 X2S + tk 2 X2 2 - (1 - x)(l - y) 

and therefore 



Xfei / 2-(l-.T)(l-j/) _ ^ / 2-(l-x)(l + y) 

- Vi V 2 - (1 - + y) ' " V 2 - (1 - - 2/) ^ ^ ^ 



and 



dxi dx2 = dMl^ dYzz. (2.25) 

Xo 

The range for the variable x is restricted by the requirement that both xi and X2 be less 
than 1; this gives 

^min <x<l, (2.26) 



with 



Xmin — max 



2(i + y)x2, 



(1 + X2 )2(1 _ y)2 + 16 2 + (1 _ _ ^2 



2(l-j/)x2 



62 



y/(l+x2,)2(l+y)2_16yx2, + (1 + t/)(l - xg,) ^ 



(2.27) 



Note that Xmm depends explicitly on y, and implicitly on MJ^ and Yz^ through a;(,i,X(,2- 
It can be checked that Xmin is always larger that M^^/S, as required by the definition of x. 

In the center-of-mass frame of the ZZ system, the four-momenta of the produced Z 
bosons can be parametrized in terms of two angles 6i,92- 

ki = — ^ (1, /3 sin ^2 sin^i, /3cos^2 sin^i, /3cos ^i) 

k2 = (l,-/?sin6'2sin6'i,-/3cos6'2sin6'i,-^cos6'i), (2.28) 
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with (3 given in eq. (2.8). Both Oi and 62 range between and vr. Thus, in addition to the 
Born variables v, we have now the three radiation variables r = ^2}, with 

a;min<a;<l, -l<y<l, 0<6'2<7r. (2.29) 

Following ref. [5] we define the corresponding integration measure 

d^r = dxdyde2. (2.30) 

Prom the computation of ref. [9] we obtain 

da'^^^ = o?$„ d^r X] i^ggi^^ A*) + ^qgi'"' A*) + Rgqiv, r, fi)] , (2.31) 



where 

Rqq{v, r, 11) 

Rqg{v,r, n) 
Rgqiv, r, m) 



1 



/5 



1 



1-y 



fq^ (Xl , H) /f ^ {X2 , fi) fqq{x, y, 6'l , 6*2 , M) 
1 (3 /I 



+ 



+ 



1 + y 



(2.32) 



1 l3 ( \ 



(47r)2 327r2 Ml^S \\-y) 



fq^ ' {X2 , m) /qg {x , y , 01, 02, fl) 

fo^ {xulA fq''{x2,lj) fqg{x, -y, 01, 02, n). 



The functions Rqq, Rqg, Rgq denote the regularized real emission cross sections for the 
different subprocesses. The functions fqq and fqg are regular in the limits of soft (x = 1) or 
collinear {y = ±1) emission; they arc given explicitly in eqs. (2.26,2.66) and Appendix C 
of ref. [9]. The distributions 1/(1 — x)p and 1/(1 ± y)+ are defined by 



(2.33) 
(2.34) 



As shown in ref. [9], the remnants of the collinear subtraction must also be added to 
get the full cross section. This contribution has the form 



da^""^ = d<Py dxdy^^ {L'^qi^, x, fj,) + Lgq{v, x, /x)) d{l - y) 



+ {Lqqiv, X, ll) + Lqg{v, X, /J,)) 6{1 + y) 



(2.35) 



where 



Ltg{v,x,n) 



2tt IQttS 



log 



Ml 



xn^ \1 — X 



+ 2 



log(l — x) 
1-x 



(l + x2) + (l-x) 
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X M^^^ (ML, t) f^^ {x,,/x, fx) f^' (x,2, fi) 



H21 



(2.36) 



L Jv, X, fi) 



Lgq{v,x,n) 



Lgg{v,X,fi) 



2tt IQttS 



X Mf^ {Ml„t) /f ^ (xw, ^) {x,2/x, 11) 



log 



ML 

XfJ,'^ 



1-x 



+ 2 



log(l - x) 
1-x 



{l + x^) + {l-x) 



log- 



M2 



27r 167r5 

X (ML, t) /f 1 (a;6i/x, /x) /f ^ (x,2, /x) 



+ 21og(l - x) (x^ + (1 - a;)2) + 2x{l - x) 



QQ 

27r WttS 
/(b). 



log- 



M2 



,T/7/ 



+ 21og(l - x) (a;2 + (1 - a;)^) + 2a;(l - x) 



X (ML, t) /f (X6i, /x) /,^^ (x;.2/x, /x) 



(2.37) 



(2.38) 



(2.39) 



and t is given in eq. (2.7). From cq. (2.27) we see that the integration range becomes 
Xbi<x<l for eqs. (2.36,2.38) and Xb2 < x < 1 for eqs. (2.37,2.39). 



3. Hcirdest event cross section 

In this section, we will write eq. (1.2) for the case at hand in full detail. The method 
of ref. [5], when applied to a generic process, may require a separated treatment of each 
singular region. In this case this is not needed. Our choice of variables v, r is adequate for 
both collinear regions at the same time, the only difference being the sign of y. We have 
instead to pay attention to the flavour structure of the process. In ordinary SMC codes, 
the flavour structure of the Born subprocess is not altered by subsequent radiation. On 
the other hand, if the hardest radiation is produced in the context of a NLO calculation, 
the association of the NLO process to a Born subprocess is not always obvious. A given 
real-emission contribution must be associated to the Born process in which it factorizes in 
the collinear limit. In the present case, the collinear regions for the qq subprocess always 
factorize in terms of the qq Born process, and the same holds for the collinear regions of the 
gq and qg subprocesses, as shown in fig. 1. Thus, for a given flavour q, we lump together 
the qq, the qg and the gq real-emission subprocesses. 

Following ref. [5], we write the cross section for the hardest event as 



da = ^ Bg{v, Hv)d^v 



(3.1) 

where the on top of R means that we drop the -|- prescriptions that regularize the x and y 
singularities. The R are thus the unregularized real emission cross sections (corresponding 
to R in the notation of ref. [5]). Furthermore, 



Bg 



{V,H) = Bq{v,n) +Vq{v,fl) + J d^r [Rqq{v,r,n) + Rqg{v , V, fx) +Rgq{v,r,n)] 
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Figure 1: NLO subprocesses that share the same elementary flavour structure. 




1 



+ 



dx [L+g{v,x,fi) + Lgq{v,x,n)] S{l~y) 



+ 



dx [Lgg{v,x,fi) + Lqg{v,x,n)] 6{l + y) 



(3.2) 



^^(pt) = exp - 



/ 



Rqqiv, r, fir) + Rqgjv, r, fir) + Rgqjv, r, Mr) 
Bq{v,Hr) 



6{kT:{v,r) -pT)d^r 



(3.3) 



and kT{v,r) is the transverse momentum of the radiated parton, 



kT{v,r) 




(l-xni-y^). 



(3.4) 



Equation (3.1) is the analogue of eq. (5.10) of ref. [5]. The function Aq{pT.) corresponds to 
A^^^^\pt) in the notation of ref. [5]. 

4. Scale choices and Sudakov form factor 

The scale choices in eqs. (3.1-3.3) deserve particular attention. We denote by a scale 
that depends only upon the Born variables, and by fir a scale that is appropriate to 
the radiation process, that is to say, of order kj-- Appropriate choices are, for example, 

fly = Mzz and fir = k^{v,r). Thus, the distinction between these two scales becomes 
particularly important when kT{v,r) <C Mzz, which is in fact the region where most 
radiation is produced. The correct scale choice is important here to maintain leading log 
(LL) accuracy in the Sudakov form factor, eq. (3.3). 

We now remind the reader how the log counting is done in the Sudakov exponent. We 
call L the large logarithm log((5^/p^), where Q is a scale of the order of the upper limit for 
k-Y- The dominant terms in the exponent have the structure L{as{Q'^)L)^ , k = I, . . . ,oo. 
We call these the LL terms. The NLL terms are of order [agL)^, the NNLL terms as{asL)^, 
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and so on. Roughly speaking, in the Sudakov exponent, the dx/{l — x) singularity and the 
dy/{l — y'^) singularity both contribute a factor of L. Expanding asik-r) as 

-)2 \ ^' 



asik^) = as{Q') 



l + E(^oa.(Q')log^J 



(4.1) 



we see that LL terms are generated when both the x and y singTilarities are present, and 
NLL terms are generated when only one of them is present. Terms with no singularities 
are of NNLL importance. 

We now show how, by a suitable slight modification of a^, we can achieve next-to- 
leading (NLL) accuracy in the Sudakov form factor. To see this, we isolate the singular 
part in the integrand of the Sudakov exponent 

Rqq{v,r,IJ,r) + Rqg{v,r,ljLr) + Rgq{v,r,ljLr) _ Cyas I 1 + x'^ f^^{Xbl/x,llr) 



+ 



Bq{v,IIr) 

T^as 1 x2 + (1 - x2) 1 (Xw/X, ^Xr) 



27r TT 



y 



+ 



1^2 



27r 7r(l-y)(l-x) x f^\x,^, iXr) 

+ regular terms . (4.2) 



The above equation is a consequence of coUinear and soft factorization, and can be easily 
obtained from eqs. (2.26), (2.42), (2.61) and (2.66) of ref. [9], together with the definitions 
of B and i?, eqs. (2.13) and (2.32) in this paper. We now replace 



fq''{Xbl/x,^r) 
xf^^{Xbl,l^r) 



fq^{xhl/x,^r) 



Xfq^ {xi,i,Hr) 



+ 1 



(4.3) 



so that the term in square bracket vanishes for x 
variable, trading y for fex, according to the formula 



1. Then we perform a change of 



1.2 

/Vrr-i 



l-J/2 



4:xk^ 



(1 



x)2M|, 



(4.4) 



and work out the Sudakov exponent up to NLL accuracy. We obtain 



lOgAg(pT) ^ 



-J 

J p. 





CpOsikl) 


U2 

A/FJl 


27r 




nasikl) 


1.2 


27r 




Cj^asik^) 



, l + x^ 

dx 

'^-x 



f^Hxbi/x,kT) 



xfq ^ {Xbl , ki 



+ 



/ dx [, 
Jo 



TT 



2 , /I ^2l f^^i^bi/x,kT) 
x^ + [1- xY\ -^-^ 

xfq ^{xbi,kT) 

dx 1 + 



+ 



-I 



dkldlogf^^{xbi,kT 



d log k^ 

dkl Cpa,(fc2) 



dHdlog/f2(xb2,A:x) 



TT 



d log 



1^2 

y^-y 

1^2 
q^q 

y*^-y 



(4.5) 
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The above formula has been obtained with the following assumptions 

• The scale Q is a scale of the order of the upper limit for /c^- Its precise value affects 
the result only by terms of order a^, i.e. next-to-next-to leading logarithmic (NNLL) 
terms. 



We assume an implicit theta function associated to the pdf 's, that sets them to zero 
when the parton momentum fraction is greater than 1. 

We replace 

dk^ 2ydy dy 



kl ■ i±y (^-^^ 

in the terms with no singularities in the x integration, and y singularities in the 
y = ±1 regions, the error for the replacement being of NNLL order. 

The upper limit of the x integration can be set to 1 in the integrals that do not have 
an X — ^ 1 singularity. In the terms singular for x 1, we replace 



V (l-x)2M|, V (l-x)2M|,' ^^-'^ 

and thus set the upper limit of integration in x to 1 — ik-r/Mzz-, since this change 
makes only subleading differences. 

• We have used the leading order Altarelli-Parisi equations. Subleading corrections to 
the evolution yield NNLL correction to the exponent. 

We obtain 



. f^'ixbi,PT)ft(^b2,PT) j f"^" dk^C^asikl) 



I 



(4.8) 



Equation (4.8) corresponds to the NLL expression of the Sudakov form factor in the DDT 
formulation [12]. In fact, Ag multiplies the Born term, that includes parton density func- 
tions evaluated at a scale of the order of Q. The double ratio of parton density functions in 
eq. (4.8) thus replaces the parton density functions included in the Born term with those 
evaluated at the scale p^, as required in the DDT formulation. The exponent in eq. (4.8) 
corresponds to the DDT exponent. It is easy to check that, by replacing 

as{kl) ^ as{kl) + i- - TT^ - al{kl) (4.9) 

in eq. (4.8) we automatically achieve NLL accuracy in our Sudakov form factor^ [13, 14]. 
In fact, we may as well perform the replacement eq. (4.9) in our initial expression for R in 
eq. (3.3). The remaining effects of such replacements in the derivation carried above are 
in fact at the NNLL level. 

The replacement eq. (4.9) was also advocated in ref. [15], in the equivalent form of an 
effective A to be use in shower Monte Carlo programs, Amc = 1-569Ams, in the framework 
of the resummation of large log(l — x) effects in the threshold region. 

^In the language commonly used in Sudakov resummations, this amounts to the inclusion of the A2 
term, that arises from the most singular contribution of the next-to-leading order Pqq splitting function. 
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5. Hardest event generation 



A straightforward implementation of the hardest event generation based upon the hardest 
event cross section cq. (3.1) and standard Monte Carlo techniques would not be very 
practical. In fact, the prefactor Bq in eq. (3.1), as well as the Sudakov exponent, already 
require a three dimensional integration over the radiation variables. This may not be a 
severe problem in the case of Z pair production, since the production formulae are relatively 
simple, but may become prohibitive for more complex processes. In the following, we will 
illustrate a method for generating hardest events with high efficiency. This method is quite 
non-trivial, and it demonstrates the applicability of the approach of ref. [5] in the most 
problematic case when initial state radiation is present. 

For ease of presentation, we first illustrate how the generation of hard events according 
to eq. (3.1) would proceed with standard Monte Carlo techniques. The sequence is as 
follows: 

1. The total cross section atot is computed as 

o-tot = X] / d^vBq{v,fiy) . (5.1) 

2. Random values of v and of the flavour label q are generated with probability distri- 
bution Bq{v,fXv), using standard hit-and-miss techniques, schematically described in 
Appendix A (see also ref. [16]). 

3. The radiation variables are generated as follows: a real random number < n < 1 is 
generated uniformly, and the equation 

n = i^qip^) (5.2) 

is solved for p-^. If there is no solution (i.e., if the resulting j»x is below the infrared 
threshold) no radiation is generated, and the event is produced as is. Otherwise, 
radiation variables r are generated with a distribution proportional to 

where p (for "process") stands for qq, qg, or gq. To be more specific one can use the 
6 function to eliminate one variable, for example x, compute 



Dp{v, y, 62) = / dx 6{k^ - p^) ^ = 

J Bq{v,Hr) 2^Bq{v,fJ,r) 



(5.4) 



where x is such that kT{x, y, v) = Pt, and then generate y, 62, and p values distributed 
with a probability proportional to Dp{v,y,62) with hit-and-miss techniques. 

The events generated according to the above prescription have uniform weights, given by 
the total cross section computed at step 1 divided by the total number of generated events. 
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The procedure outlined above is computationally intensive, both in the generation of 
the Born variables, which requires one or more three-dimensional integrations per generated 
point, and in the generation of the radiation variables, where the computation of the 
Sudakov exponent also requires a three-dimensional integration. We shall now illustrate 
our method. We shall discuss separately the generation of Born variables and the generation 
of radiation variables. 

5.1 Generation of the Born variables 

We first replace the radiation variable x by a rescaled variable x ranging between and 1: 

* ^ — a^min _ _ x fr. r,\ 

X ^ , X — S^min -p XyL 3^minj) V"""^/ 

where Xmm = Xmmiy, v) is given by eq. (2.27). In the two collinear regions y = ±1 we have, 
respectively, Xrain{^,v) = Xbi and Xmin(— l,^') = Xb2. Therefore, we also define 

X~^ = Xbl + X{1 - Xbl) , X~ = Xb2 + X{1 — Xb2) ■ (5.6) 

The new radiation variables r = {x,?/, ^2} have now fixed integration ranges 

0<x<l, -l<y<l, 0<6'2<7r. (5.7) 

We rewrite Bq{v,iJ,) as 

+ J dy (162 dx {1 - Xmin) [Rqq{v,r,yL) + Rqg{v,r,iJL) + Rgq{v,r,ix)\ 

+ J dx{l-Xbl)[L'^qiv,X^,fJ,)+Lgq{v,X^,fJ,)] 

+ J dx{l-Xb2) [L-g{v,x~,n) + Lqg{v,x-,n)] 

dy d62dx Bq{v,f, n) , (5-8) 



where 



Bq{v,r,n) = ^[Bq{v,lj) +Vq{v,lj)] 



+(1 ) [Rqq{v, r, IJ) + Rqg{v, r, jj) + Rgq{v, r, //)] 

1 

1 
'2^ 



+^{1-Xbi) [L-^g{v,x+,n) + Lgq{v,x+,n)] 
+-^{^-Xb2) [Lgg{v,x',n) + Lqg{v,x-,ii)] , (5.9) 



and we define 



B{v,r,ti) = J2Bq{v,f,n) , (5.10) 
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so that 



atot = J d<^rd<PvB{v,f,iJ,y), (5.11) 



where d^f. = dy d92 dx. 

We now observe that the integrand in eq. (5.11) is, in general, positive. It is in fact 
the sum of a positive term of order in as and terms of order ag- Negative terms are 
not totally forbidden, but they clearly disappear in the perturbative (i.e. ag — ^ 0) limit. 
It will be thus necessary, when performing the numerical integration, to check that the 
contribution of negative terms is negligible. 

It is possible to store certain intermediate results of the numerical integration pro- 
cedure, that can be subsequently used to generate efficiently the variables v,f with a 
distribution proportional to B{v,r,^y). We give in appendix B an elementary description 
of such a procedure. In practice, computer programs are available that implement and 
optimize this procedure. We have used the BASES/SPRING [17] package. The adaptive 
Monte Carlo integration routine BASES performs the integration and stores the necessary 
intermediate results. The routine SPRING uses this information to generate unweighted 
events according to the distribution B (v , f , fj,^) . In this way, both v and f variables are 
generated, and then the flavour q of the event is generated with a probability proportional 
to Bq{v,r,Hy). The values of the radiation variables r are then ignored, since we are only 
interested in the distribution of the Born variables v. This precisely amounts to integrat- 
ing away the r variables, so that we are left with a uniform generation of v and q values 
according to the Bq(v,fj,) distribution. 

5.2 Generation of the radiation variables 

Having generated the Born variables v and the flavour q, the next step is the generation of 
the radiation variables according to the probability distribution 

Aq{v,k^{v,r))Wq{v,r)d^r , (5.12) 

where 

^(^^ ^ ^-/^('■■''-/'' ) + ^^^(;'''-/'; ) + ^^^('''^'/^-) (5.13) 

Bq{v,IJ,r) 



\iPT,v) = exp 



J Wq{v,r')0{k.r{v,r')-p^)d^r 



(5.14) 



We do this by means of the veto technique, described in some detail in appendix C. We 
find a function Uq{v,r) such that 

Uq{v,r)>Wq{v,r) , (5.15) 
and we generate radiation variables r according to the distribution 

A'f\v,kr{v,r))Uq{v,r)d^r (5.16) 

where 



J Uq{v,r')9{k^{v,r')-p^)d^r 



(5.17) 



The generation of the event is then performed by the following steps: 
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1. Set Pmax = OO. 



2. Generate a random number n between and 1, and solve the equation 



n= \ (5-18) 



for 

3. Generate the variables r according to the distribution 

Uq{v,r)6ik^{v,r)-p^) . (5.19) 

4. Accept the generated value with a probability Wq{v,r)/Uq{v,r). If the event is re- 
jected, set Pmax = Pt, and go to step 2. 

It is obvious that an efficient generation is achieved if the chosen functional form for Uq is 
such that the solution of eq. (5.18) and the generation of r variables at step 3 are simple, 
and that the ratio Wq{v, r)/Uq(v, r) is not too far below 1 in most of the integration range. 
We find that the choice 

U(vr)-N (5 20) 

with a suitable choice of the constant Nq, fulfills both requirements. The generation of 
events according to the distribution in eq. (5.16), with Uq(v,r) given by eq. (5.20), is not 
entirely trivial; we describe it in Appendix D. 

The choice in eq. (5.20) is suggested by the structure of the function Wq in eq. (5.13) 
near the collinear limit, as we now show. Let us consider first the qq contribution. In the 
positive collinear direction (y — 1) we have (see eq. (4.2)) 

Rqqiv,r,l2r) ^ 1 1 CpQ^ 1 + f^^{.Xbl/x,fXr) 
Bq{v,Hr) ^TTl-y 27r 1-X xf^'{Xbl,l^r) ' 

A similar relation holds for y ^ —1. It is now reasonable to assume that the parton 
density ratio in eq. (5.21) is bounded by a number of order one, since parton densities, in 
general, are never fast growing functions of x. A bounding function of the form of eq. (5.20) 
therefore arises in a natural way. We now consider the gq contribution. In the y — ^ 1 limit 
we have 

^ [(1 - x)^ + X-] f^'^^^^l-^^'^^ . (5.22) 

In this case, the ratio of parton densities is between different parton species, and must be 
discussed with care. First of all, we notice that 

/f 1 (Xbi/X, Mr) < fg' (xw, Mr) , (5.23) 
so that we only need a bound for 

ll^^^. (5.24) 

fq \Xbi,IJ,r) 
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At small values of Xhi, this ratio is always bounded, because the gluon and quark densities 
have similar behaviour in the small-x limit. For large x^i, if g is a valence quark, the ratio 
is also bounded, since the gluon is generated by valence quarks through evolution. In case 
g is a sea quark, the corresponding density may be softer than the gluon density. In the 
worst case, however, the sea quark is generated by the gluon through evolution. Assuming 
that parton densities behave as a power of 1 — x at large x, 

fg{x,n)^{l-x)^ (5.25) 

the Altarelli-Parisi equation in the large x limit yields 

ajj,'^ 27r z 27r 



and therefore 



Since 



f^^{xbi,l^r) / 1 



< 



(5.27) 



11 , . 

< , 5.28 

1 — Xfti 1 — a:; 

we conclude that the choice eq. (5.20) is adequate also in this case. 

The normalization factor Nq can be obtained numerically, using SPRING to generate 
u, r and q values, and then computing the maximum of the ratio Uq{v, r)/Wq{v, r). 

5.3 Colour assignement 

The NLO calculation of a generic production process may specify also the detailed colour 
structure of the produced particles. SMC generators use the colour information for the 
subsequent shower, and also at the hadronization stage, for the formation of colour singlet 
hadrons. In general, only the colour flow structure in the limit of a large A^c (where Nc is 
the number of colours) is needed, and, in fact, in the Les Houches Interface one can only 
specify the colour connections of the intervening partons. In the case of ZZ production, 
the large colour assignment is the same in all real emission graphs, corresponding to 
the quark and antiquark colour line merging into the gluon double line. This is thus the 
colour structure that must be passed to the LHI. 

In the general case, several colour configurations are possible, and one should specify 
which one to choose after the radiation has been generated. If the contribution to the real 
emission cross section is available for each (large Nc) colour component, one simply chooses 
the colour component with a probability proportional to it. 



6. Results and comparisons 

We now illustrate some numerical results obtained by the hardest event generation prescrip- 
tion presented in the previous Sections. We will refer to this method as as the POWHEG 
(for Positive weight Hardest event Generator) . Our results are obtained with the following 
default choices: 
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• We use the CTEQ6M [18] parton density set. 

• We fix the factorization and renormalization scales in the computation of Bq to the 
invariant mass of the Z boson pair, fj,^ = Mzz- 

• We fix the factorization and renormahzation scales for radiation to the k-j- of the 
radiated parton, eq. (3.4). 

• We use the value of A^s appropriate to the PDF set we have chosen. This value is 
corrected according to the prescription given in Section 3 when generating radiation 
variables. 

• We consider two experimental configurations of interest: pp collisions at s/S = 
1920 GeV, corresponding to the Tevatron, and pp collisions at s/S = 14 TeV, corre- 
sponding to the LHC. 

• We use a fixed number of flavours nf = 5. In principle this choice is not completely 
consistent. One should instead reduce nf when generating radiation below the bottom 
and charm scales. This has however a minor impact on the final results, and we have 
chosen not to take it into account at this stage. 

• The Z bosons are treated as stable particles, with = 91.2 GeV. We have forced, 
whenever possible, the same assumption on standard Monte Carlo predictions, when 
comparing them to our results. 

• We have used the following values for the electroweak parameters: sin^ Ow = 0.23113 
and aem = aem{ml) ~ 1/128. 

We begin by comparing the POWHEG results with the fixed order QCD calculation 
[9], performed with the same scale choice adopted for Bq, fi = Mzz- We have considered 
distributions in the following variables: the transverse momentum and rapidity of one Z 
boson, the rapidity and pseudorapidity differences Ay and Ar/ between the two Z's, the pair 
invariant mass Mzz, the transverse momentum and rapidity of the ZZ pair, the azimuthal 
distance A0 between the two Z's, and Ai? = ^/A^~+~Arp. Among the distributions 
we have considered, the A0 and p^z distributions are strongly affected by light parton 
emission, since they are trivial in the case of no emission. The AR distribution is an 
intermediate case, since it depends upon A<p. All other variables are non-trivial already at 
the Born level. We find that the NLO calculation and POWHEG give equivalent results 
for this last group of observables, while in the case of A(f) and pJz' and, to a lesser extent, of 
AR, we find important differences. This is illustrated in fig. 2, where the A^, p^^ and AR 
distributions are shown, together with the invariant mass distribution, which is taken as a 
representative example of all the remaining observables. We first notice that in the case of 
the invariant mass distributions the two calculations give identical results. On the contrary, 
the A(/) and p^z distributions are strongly affected by light parton emission, and indeed 
they display a sizable difference in the two calculations. In the fixed order calculation the 
divergences in the limit of soft/collinear light parton emissions are canceled when virtual 
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Figure 2: Comparison of four distributions computed according to the fixed order calculation and 
the POWHEG. On the vertical axes we report the cross sections in picobarns per bin. Energies are 
expressed in GeV. 



corrections, that have strictly p^z = and Ac/) = vr, are included. However, the effect of 
virtual corrections is not seen in the At/) and p^z plots, so that the fixed order calculation 
appears to yield logarithmically divergent integrals for p^z ~^ ^'^d A^ — > vr. On the other 
hand, the POWHEG result is well behaved also in this region, since the no-radiation region 
is effectively suppressed by the form factor Aq{kT^) (see eq. (3.1)). 

Our next task is to compare our full result with the only existing implementation of 
Monte Carlo generation of vector boson pairs with NLO accuracy, namely the MC@NLO 
program. The relevant plots are collected in figs. 3, 4, 5 for the Tevatron, and in figs. 
6, 7, 8 for the LHC. The MC@NLO results are compared with POWHEG interfaced with 
HERWIG. According to the Les Houches interface prescription, the showering in the Monte 
Carlo is vetoed by assigning the /ct of the event generated by POWHEG to the variable 
SCALUP in the Les Houches common block. We find that the inclusion of the HERWIG 
shower determines only tiny changes with respect to the pure POWHEG output. It is easy 
to comment upon the outcome of this comparison: the two algorithms yield identical results. 
Minor differences are only seen in the p^z distribution; these can be easily attributed to 
the presence of a k-r hard cutoff, set at 0.8 GeV in the POWHEG. 



-19- 



It should be noted that the MC@NLO results have been computed using the scale 
choice suggested by the authors, namely 

/^^ = ^ + + \/pt2 + "^|) ' (6-1) 

while the POWHEG results are obtained with the scale choices described at the beginning 
of this section. This difference has a negligible impact on the total cross section, but 
may affect some of the distributions. The agreement between the two methods in the 
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Figure 3: Comparison of four distributions computed according to MC@NLO and the POWHEG. 
On the vertical axes we report the cross sections in picobarns per bin. Energies are expressed in 
GeV. 

region of large transverse momentum radiation is expected, since both methods are in good 
agreement with the fixed-order calculation in this region. The region of soft radiation in the 
MC@NLO implementation is controlled by the HERWIG shower. We therefore conclude 
that our treatment of the soft region is consistent with the HERWIG implementation. We 
point out that HERWIG fully implements the procedure discussed at the end of Section 3 
for the treatment of the strong coupling constant in the shower emission, as we also do. 

One of the main features of the POWHEG method is the possibility of interfacing its 
output to any shower Monte Carlo that implements the Les Houches interface for user- 
provided processes. This is possible with the popular Monte Carlo PYTHIA [19] since 
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Figure 4: Comparison of four distributions computed according to MC@NLO and the POWHEG. 
On the vertical axes we report the cross sections in picobarns per bin. Energies are expressed in 
GeV. 

its version 6.3, which is based upon a px-ordered shower algorithm. We thus interfaced 
POWHEG to PYTHIA and compared the results to the pure PYTHIA output, normalized 
to the POWHEG cross section. By inspection of fig. 9, relevant to the Tevatron configura- 
tion, we see that the POWHEG result appears to have a harder p^z spectrum, especially 
in the low p^z region. Correspondingly, the Ac/; distribution obtained by POWHEG is also 
shifted away from the region A(/> ~ tt with respect to the PYTHIA result. A detailed 
view of the small-p2z region is shown in fig. 10 for both the Tevatron and the LHC. We 
see that the position of the peak in the p^z distribution is rather different in the two ap- 
proaches. We have verified that also in this case, the effect of the PYTHIA showering in the 
POWHEG-I-PYTHIA result is negligible in the distributions we are considering. We thus 
infer that the same differences should appear when comparing PYTHIA with MC@NLO. 
This issue deserves further investigations, that are however beyond the purposes of the 
present work. 

7. Conclusions 

We have presented an explicit implementation of the method presented in ref. [5] for the 
construction of a Monte Carlo event generator with matrix elements accurate to next-to- 
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Figure 5: Comparison of four distributions computed according to MC@NLO and the POWHEG. 
On the vertical axes we report the cross sections in picobarns per bin. Energies are expressed in 
GeV. 



leading order, and positive weights. The process we have considered is Z pair production in 
hadron colhsions, which is a process of intermediate complexity that involves initial-state 
radiation, and thus is a good testing ground for the applicability of the method to hadron 
collision processes. 

One of the main achievements of the present work is the development of numerical 
techniques that are necessary for the practical implementation of the method proposed in 
ref. [5]. We have given here a full, detailed description of these techniques, so that all our 
results can in principle be reproduced by the interested reader. 

The output of our generator is cast in the form of the Les Houches Interface for user 
provided processes [8], and thus easily interfaced to both the HERWIG and the PYTHIA 
shower Monte Carlo programs. 

We have compared our result to the only existing program that can compute NLO 
corrected Shower Monte Carlo output in hadronic collisions^, for a large set of observables, 
both for the Tevatron and the LHC. We have found an excellent agreement. 

Our method can also be used as a standalone, alternative implementation of QCD 

* Other methods for NLO implementations for shower Monte Carlo have been presented in the literature 
[20], that are limited to the case of e+e~ annihilation, and do not have positive weights. 
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Figure 6: Same as fig. 3 for the LHC. 



corrections. As discussed in section 6, the interfacing to the SMC has neghgible effects on 
the POWHEG distributions involving ZZ observables. It may thus becomes convenient to 
compute these observable by using POWHEG as a standalone program. This has several 
advantages over the standard, fixed order QCD calculation, like positivity, and next-to- 
leading logarithmic accuracy in the small /ct region. 

The extension of the present work to W Z and WW pair production is straightforward. 
Furthermore, the choice of kinematic variables adopted here can be easily extended to 
processes of Drell-Yan type, like single W and Z production, and Higgs production. We 
also believe that the method can be applied without major modifications to the case of 
heavy flavour production, where final state soft radiation is also present. The treatment of 
collinear final state radiation with the method of ref. [5] presents no difficulties, being in 
fact easier than the initial state radiation case. The formulation of the application of the 
method to general processes is under study. 

In order to fully preserve double logarithmic accuracy in the POWHEG context, the 
Shower Monte Carlo to which the program is interfaced should generate consistently the 
needed soft radiation pattern. When using an angular ordered shower program, one needs 
to include additional soft radiation, that was named "truncated showers" in ref. [5]. Pro- 
grams that generate soft radiation using dipole type formalisms [6, 7] may not require any 
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Figure 7: Same as fig. 4 for the LHC. 

further correction. At present, however, no such program is available that handles initial 
state radiation. These problems are left for future studies. 

Acknowledgments: We thank Stefano Frixione for helpful discussions and suggestions, 
and Carlo Oleari for carefully reading the manuscript. We also thank Torbjorn Sjostrand 
for clarifications on the use of PYTHIA. 

A. The hit-and-miss technique 

In order to generate a set of continuous variables x and an index j = 1, . . . ,m, with a 
probability distribution proportional to Fj{x), one first computes 

Ftot = ^ dx Fj (x) , Fmax = max ^ Fj (x) . (A.l) 

j j 

Then, one generates a random value x = x and a random number < n < 1 with flat 
distributions, and finds the smallest value k such that 

k 

Y,Fj{x)>nF,^^^. (A.2) 
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Figure 8: Same as fig. 5 for the LHC. 
If such k does not exist, the event is rejected; otherwise the set j = k, x = x is accepted. 

B. Prom Monte Carlo integration to uniform generation 

We want to to generate a set of values for the variables x in a domain D with a distribution 
proportional to F{x). We divide the integration range P into a set of hypercubes I?,, with 
i = 1, . . . , m. For each hypercube we compute 

„ m 

fI^,= dxF{x), F41 = maxF(x), F,,, = Y,Fi^- (B-l) 

In order to generate the x values, we first choose the hypercube according to its probability: 
given a random number < n < 1 we find the minimum value k such that 

k 

5^Ft(^^>nFtot. (B.2) 

i=i 

Next we generate x in the Dk hypercube using the hit-and-miss technique. It is clear 
that the efficiency of the generation improves by reducing the size (and thus increasing 
the number) of the hypercubes. Thus there will be a trade-off between speed and storage 
requirement in the computer implementation of this technique. 
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Figure 9: Comparison of four distributions computed according to PYTHIA and the POWHEG 
interfaced to PYTHIA. 
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Figure 10: Comparison of the p^z distribution computed according to PYTHIA and the POWHEG 
interfaced to PYTHIA, for both the Tevatron and the LHC. 



C. The veto technique 

Assume we want to generate values for a set of variables x, according to a distribution 



P{x) = R{x) exp 



d'^x' R{x')9{p{x') -p{x)) 



(C.l) 
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where p{x) > 0. We assume that R{x) is non-negative, and that the unconstrained integral 
is divergent. We observe that the probabihty distribution of p{x) is an exact differential. 
Indeed, 

J d'^x P{x) S{p{x) -P) = J d'^x R{x) 5{p{x) - p) exp - j R{x') e{p{x') - p) 



where we have defined 



A(p) = exp 



/A'i,MWx')-p) 



(C.3) 



The function A(p) ranges between and 1, because the integral in the exponent vanishes 
for p +00 and diverges to +00 for p = 0. Hence, eq. (C.2) also shows that P{x) is 
normahzed to 1: 

//•+00 /• 
d'^x P{x) = J dp j d'^x P{x) 5(p{x) -p) = A(oo) - A(0) = 1 . (C.4) 

In principle, the uniform generation of events is therefore straightforward: one generates 
a uniform random number r between and 1, solves the equation A(p) = r for p, and 
then generates values of x on the surface 5{p{x) — p) with a distribution proportional to 
R{x)5{p{x)-p). 

In practice, however, the solution of the equation A(p) = r is in most cases very heavy 
from the numerical point of view. This difficulty can be overcome by means of the so-called 
vetoing method, which we now describe. We assume that there is a function H{x) > R{x) 
for all X values, and that 



Aif(p) = exp 



-/A'i^M.«x')-p) 



(C.5) 



has a simple form, so that the solution of the equation Ah{p) = r and the generation of 
the distribution H{x) 5{p{x) —p) are reasonably simple. Then, we implement the following 
procedure: 

1. Set Pmax = 00. 

2. Generate a flat random number < n < 1. 

3. Solve the equation 

^^^^^=n (C.6) 



Afl-(pmax! 

for p (a solution with Q <p < Pmax always exists for < n < 1). 

4. Generate x according to H{x) 5{p{x) — p). 

5. Generate a new random number n' . 
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6. If n' > R{x)/H{x) then the event is vetoed, we set Pmax = go to step 3 and 
continue; otherwise the event is accepted, and the procedure stops. 

The resulting events are distributed according to eq. (C.l). The proof of this statement 
is simple but non-trivial. At the end of the vetoing procedure, the event distribution will 
be given by the sum of the distribution for the case in which there is no veto, there is one 
veto applied, two vetoes, etc.. The probability distribution of events generated with no 
veto applied is given by 



JO 



dpi . \ H{x) 5{p{x) - pi)) = R{x) Ah{p{x)) . 



^niPn 



H{x) 



(C.7) 



We have used the fact that AHipma^) = Ah{oo) = 1, and we have inserted a factor of 
R{x)/H{x), corresponding to the acceptance probability. 
When one veto is applied, we have 



Pi(x) 



dp 



Ah{pi 



Ah{p, 

Ah{P2) 

Ah[pi) 



d'^xiH{xi)5ipixi)-pi) 1 



Rixi) \ 
H{x,)J 



VI 



H{x)5{p{x)-p2)) 



R{x) 
H{xj 



+00 



R{x)AH{pix)) / dpih{pi), 

Jp{x) 



(C.8) 



where we have defined 



hipi) = J d'^xiH{xi)5{pixi) -pi) 



Hixi)) 



(C.9) 



The factor 1 — R{x\)/H{x\) is the rejection probability, which must be inserted at each 
vetoed step. Note that the result is nonzero only for pi > p{x), because of the 6 function 
in the p2 integration. It will be useful to perform one more step explicitly; for two vetoes, 
we find 



Jo 



dpi 



Ah{p, 

Jo ^ Ah{P2) 



■h{Pi) / dp2 . , . h{p2) 

Jo ^H[Pl) 



H{x)S{p{x)-ps)) 



R{x) 
H{x) 



R{x) Ah{p{x))]^ 



+ 00 



p{x) 



dp h{p) 



(C.IO) 



where we have used symmetric integration. It is now easy to obtain the generic term of 
this infinite sum, namely the term with n vetoes applied. We get 



Pn{x) = AH{p{x))R{x)^^ 



p{x) 



dp h{p) 



(C.ll) 
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The sum over n yields 

oo 

^Pn{x) = R{x)AHip{x)) exp 
= R{x) Ah{p{x)) exp 



n=0 



dp h{p) 

p{x) 

d'^x' [H{x') - R{x')] e{p{x') - p{x)) 



R{x) exp 



which is the announced result. 



d'^x' R{x')e{p{x') - p{x)) 



(C.12) 



D. Generation of events according to A 



The veto technique described in the previous Appendix can be employed to find the solution 
of eq. (5.18), which we reproduce here: 



n 



^f\v,PT) 
^'f\v,Pmeoc) 



(D.l) 



where n is a random number between and 1, 



Pt) = exp 



and 



Uq{v,r) e{krr{v,r) - Pt) d^r 



^ (1 - x) (1 - y^) 



r = {x,y,92} (D.2) 



(D.3) 



Trading the variable y for 



kl = kl{v,r) = ^ (1 - xfil - y% \y\ = Jl - -p^T^ (D.4) 



4x 



(1 - x)^ Ml, 



we find 



Uq{v,r)e{kT:{v,r)-pT)d^r = I dx dy de2Uq{v,r) 0{kT{v,r) - p^) 

Jp J-i Jo 

= nN, r'~'4 , °-"-?) (D.5) 



where 



M2 



' Tmax 



X±= L 1 + 



i-2 h 



(D.6) 



We have supplied an extra factor of 2 to account for the fact that each value of k?^ cor- 
responds to two values of y. The integration range for x is limited by the definition 
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X = M^^/s, and by the condition s < S. The x integration can be performed with the 
variable change = -^a;+ — x + ^X- — x. We get 



Uq{v,r)e{kT:iv,r) -pT,)d<i>r 



1.2 



where 



V{kl)=TrNgasikl) log 



y/X+ - P + - p 

^JX+ - p - y/X- - p 



We now observe that 



y/X+- p + y/X- - p + 1 , + ^Iz ^ 1 1 <f 



where 



Furthermore, 



as{kl) 



^-^a%kl)<a',{k-i), 



where ag{k^) is the leading-log expression of the running coupling, 



1 



Hence 



y(fc?)<y(fc?) = ^^iogiog|^. 



The A;^ integral of V can be performed analytically: 



Tmax 



2/3o 



l0g|2l0, 



,2 log "^"19-^ 



logS 



log^M 



(D.7) 
(D.8) 

(D.9) 

(D.IO) 
(D.ll) 

(D.12) 
(D.13) 



(D.14) 



We now proceed as follows: 

1. We set Pmax = ^Tmax- 

2. We generate uniformly a random number n, < n < 1, and we solve numerically the 
equation 



n = — — ; '[pT,)=exp 



A(^)(p^ax)' 



'^Tmax dk 



2 



^ Viki) 

T 



(D.15) 



for Pt- 



3. We generate a second random number n' between and 1; if n' < V{p'^)/V{p'^) we 
accept the event, otherwise we set Pmax = Pt and we return to step 2. 
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